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ABSTRACT 

Massive stars disproportionately influence their surroundings. How they form has only started to 
become clear recently through radiation gas dynamical simulations. However, until now, no simulation 
has simultaneously included both magnetic fields and ionizing radiation. Here we present the results 
from the first radiation-magnetohydrodynamical (RMHD) simulation including ionization feedback, 
comparing an RMHD model of a 1000 M rotating cloud to earlier radiation gas dynamical models 
with the same initial density and velocity distributions. We find that despite starting with a strongly 
supercritical mass to flux ratio, the magnetic field has three effects. First, the field offers locally 
support against gravitational collapse in the accretion flow, substantially reducing the amount of 
secondary fragmentation in comparison to the gas dynamical case. Second, the field drains angular 
momentum from the collapsing gas, further increasing the amount of material available for accretion 
by the central, massive, protostar, and thus increasing its final mass by about 50% from the purely gas 
dynamical case. Third, the field is wound up by the rotation of the flow, driving a tower flow. However, 
this flow never achieves the strength seen in low-mass star formation simulations for two reasons: 
gravitational fragmentation disrupts the circular flow in the central regions where the protostars 
form, and the expanding H II regions tend to further disrupt the field geometry. Therefore, ionizing 
radiation is likely to dominate outflow dynamics in regions of massive star formation. 



1. INTRODUCTION 

Massive stars influence their surroundings through 
their ionizing radiation, winds, and supernova explosions. 
However, their formation re mains less well understood 
than that of low mass stars ( |Mac Low fc Klessen||2004[ ) . 
Massive stars accrete much taster than low m ass stars 
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at rates r eaching as high as lO -3 M0yr (e._ 
|et al.|2002 ). This is required to build up the mass ol the 
star p rior to exhausting its nuclear fuel flKeto fc Wood 
[2006| . They are observed to have companions with far 
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greater frequency tha n low mass stars ( |Lada|[2006 
jnecker fc Yorke| 2007[ ). The ionizing radiation they emit 
produces a set of observables distinctly different from low 
mass star formation. 

Magnetic fields pervade neutral atomic gas in the in- 
terstellar medium at strengths high enough to prevent 
isotropic gravitational collapse ( |Heilesp976[ ). However, 
gravitational collapse along field lines can proceed unim- 
peded, and, in some circumstances, magnetic field s can 
even accelerate collapse through magnet ic braking (Kim 



et al.||2003| |Hennebelle fc Ciardij [2009] ) . Observations 
ol low- mass star forming regions show that magnetic 
field strength remains roughly constant at values of B = 
5-10 /iG at number densities n < 10 3 cm -3 , with the 
stre ngth growing as B oc n 2 / 3 at higher number densi- 
ties ( |Crutcher|[2009l ). 

The role ol magnetic fields in low mass star formation 
has received extensive attention. They act in three ma- 

thomas.peters@ita.uni-heidelberg.de 

1 Fellow of the Baden-Wiirttemberg Stiftung 

2 Zentrum fur Astronomie der Universitat Heidelberg, Institut 
fur Theoretische Astrophysik, Albert-Ueberle-Str. 2, D-69120 
Heidelberg, Germany 

3 Department of Astrophysics, American Museum of Natural 
History, 79th Street at Central Park West, New York, New York 
10024-5192, USA 



jor ways. First, they can provide support against gravi- 
tational collapse. Second, they drain angular momentum 
directly through magnetic braking, as well as indirectly 
through magnetorotational instability driven turbulent 
viscosity. Finally, as one consequence of angular mo- 
mentum transport, they drive jets and outflows. 

The critical mass-to-flux ratio at which magnetic fields 
can no longer support a gas cloud against gravitational 
collap se is ( |Mouschovias 
l99lt 



Spitzer 1976 Mouschovias 
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= 490 g G _1 cm- 2 



(1) 



where C = 0-53 is the value for a uniform sphere 
( Strittmatter 1966). For a field of 10 /iG, this corre- 
sponds to a column density of 5 x 10 -2 g cm -2 , well be- 
low the 0.1-1 g cm" 2 characteristic of massive star form- 
ing regions (e.g. [Mueller et al. | 2002[ |McKee fc Tan|2003 



Krumholz fc M cKee 2(M§7 As a result, less attention has 
been paid to the dynamics of magnetic fields in massive 
star forming regions where it appeared unlikely that they 
would prevent star formation (see, however, Hennebelle 
et al. 2010, A& A submitted). 

However, even in models of low-mass star formation 
starting with supercritical values of (M/<£), fragmenta- 
tion is reduced compared to the pure gas dynamical case 
by a combination of magnetic pressure and tension forces 
providing additional su pport to counteract g ravitational 
collapse (Ziegler 2005] |Banerjee fc Pudritz||2006l |Price 



|fc Bate||2007| |Henne belle & Fromang||2008 JHenne belle 
Teyssier 2008; Hennebelle & Ciardi 2009 ; Lvommercon 



|et al.||2010flBurzle et al.| |2010p. 



Magnetic fields can also drain angular momentum from 
rotating clumps and cores even if they do not entirely 
prevent collapse. Magnetic braking in idealized perpen- 
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dicular and parallel cases w as computed by Mouschovias 



& Paleologou (1979 1980), who found the criterion for 
braking to be effective is that outgoing helical Alfven 
waves pass through a mass of external gas equal to that 
of the cloud (see also ~ 



Basu & Mouschovias 



nebelle & Ciardi|2009|) 



1994 
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Magnetorotational instability further drives outward 
angular momentum transport within differentially ro- 
tating disks, or indeed any differentially rotat ing struc- 
ture with angular v elocity decreasing outward (B albus fc| 
Hawley|1991l 119981). This may be less important m mas- 



sive star forming regions where gravitatio nal instability 



dominates angular momentum transport (Peters et al. 
2010a|bp . 



Outflows have been found almost universally around 
forming low mass stars. Controv ersy persists about th e 
exact mechanism driving them ( |Ferreira et al.|E006), 
with viab l e su ggestio ns including disk win ds ( |Blandford| 

12QQQD, tower flows 



Payne] [1982] |K6nigl & Pudritz 



producea by wound up toroidal field (No misaka 1998 
2002 ] IMatsumoto fc Tomisaka|2QQ4l |Machida et al.|2UU4| 



Banerjee & Pudritz 2006, [2QQ7| ) or interactions between 



the disk fiel d and pro tostellar m agnet osphere, either at 



an X-poirit ([Najita fc Shu||1994| |Shu et al.||1994[ |2007l 
Cai et al. 2008) o r m more ge neral geometries ( |Lovelace| 
et al.||l99yl "[K omanova et al.| [2009). However, all these 
mechanisms depend on the interaction of magnetic fields 
with a coherent and well-defined rotational flow in the 
accretion disk. 

In this work we use numerical simulations to examine 
more carefully the dynamics of magnetic fields during 
the formation of massive stars, focussing on the clump 
scale of 0.1 to 0.001 pc. We model the collapse of a 
magnetized, rotating clump containing 1000 M®. This 
work is complementary to the initial numerical study 
of magnetic fields in mas sive star forming regions by 



Banerjee & Pudritz (2007) in that it examines a larger 
mass region at larger scales. It is also complementary to 
recent radiation-magnetoh ydrodynamical simu lations of 
low-mass star formation by |Tomida et al.| ( |2010[ ) because 
we take the feedback by ionizing and non-ionizing radi- 
ation into account. We show that the three ways that 
magnetic fields act in low mass star formation recur in 
massive star formation, but with different implications 
and relative importance. 

In Section [2] we describe our numerical simulations, 
while in Section [3] we discuss our results. Finally, we 
summarize and conclude in Section 0J 

2. NUMERICAL METHOD AND INITIAL CONDITIONS 

We present the first three-dimensional, radiation- 
magnetohydrodynamical simulations of massive star for- 
mation, taking into account heating by both ionizing 
and non-i onizing radiation, usi ng the adaptive-mesh code 



FLASH (Fryxell et al. 2000). We propagate the ra- 



diation on the adaptive mesh with our extended ver- 
sion of the hybrid characteristics raytraci ng method ( Ri^ 
jjkhorst et at]|2006| [Peters et al. | |2010a[) We use sm£ 
particles ( |FederratTi et al.]|2010[ ) to model young stars. 
Sink particles are inserted when the Jeans length of col- 
lapsing gas can no longer be resolved on the adaptive 
mesh. They continue to accrete any high-density gas ly- 
ing within their accretion radius. We use the sink particle 
mass and accretion rate to determine the radiation feed- 



TABLE 1 

Overview of collapse simulations. 



Name 


Radiative Feedback 


Multiple Sinks 


Magnetic Fields 


Run A 


yes 


no 


no 


Run B 


yes 


yes 


no 


Run D 


no 


yes 


no 


Run E 


yes 


yes 


yes 



back with a prestellar model ( Peters et al.||2010a ). The 
radiation-magnetohydrodynamical equations are solved 
with a novel , positive-def inite, MUSCL-Hancock, Rie- 
mann solver (Waagan 2009). 

In this work we present simulations that incorporate 
magnetic fiel ds into the gas dyna mical collapse models 



presented in Peters et al. (|2010aJ). We further analyzed 

me accr etion history of 
|2010b [), as wel l as th e 



those models with a locus on t 



the stellar cluster (iPeters et al. 



resulting H II region morphologies (|Peters et al.|[2010c). 



The initial conditions for the gas and the simulation pa- 
rameters in the simulations presented here match those 
in our previous simulations. We start with a 1000 M 
molecular cloud having a constant density core with 
p = 1.27 x 10 -20 gem -3 within a radius of r = 0.5 pc, 
surrounded by an r -3 / 2 density fall-off out to r = 1.6 pc. 
The cloud rotates as a solid body with an angular ve- 
locity uj = 1.5 x 10 _14 s _1 . The initial temperature is 
T = 30 K. The highest resolution cells on our adaptive 
mesh have a size of 98 AU. Sink particles are inserted at 
a cut-off density of p cr i t = 7 x 10 -16 gem -3 and have an 
accretion radius of r s i n k = 590 AU. 

We compare the results of four different simulations 
(see Table [I]), three of whic h were discussed and ana- 
lyzed in our previous work (Peters et al. 2010a|c|b ). In 
the first simulation (run A), a dynamical temperature 
floor is introduced to suppress secondary fragmentation. 
Only one sink particle (representing a massive protostar) 
is allowed to form. In the second simulation (run B), sec- 
ondary fragmentation is allowed, and many sink particles 
form, representing a group of stars, each contributing to 
the radiative feedback. The third simulation (run D) is 
a control run in which secondary fragmentation is still 
allowed, but no radiation feedback is included. 

The new, fourth simulation (run E) is a magnetized 
version of run B, the full stellar group simulation with 
radiation feedback from all stars. Run E includes an ini- 
tially homogeneous magnetic field along the rotation axis 
of the cloud with a magnitude of 10 /iG, corresponding 
to (M/<3>) = 14(M/$) cr in the central core at the begin- 
ning of the simulation and a plasma beta (/3 p i = Pth/.Pmag 
with the thermal and magnetic pressure p t h and p ma g, re- 
spectively) of /3 p i = 3.7 . This strength was chosen to be 
consistent with an analysis of combined H I , OH, and 
CN Zeeman measurements that yielded a field strength 



B~(10//G)( — A 

v P J V 300 cm- 3 / 



0.68 



(2) 

,300 cm- 3 J v J 

(Crutcher et al. 2010, in prep; |Crutcher|[20Tot . Our 
choice of a uniform field initially neglects the increase 
of field strength with increasing density, but this occurs 
quickly as the collaps e proceeds. We n ote that CN Zee- 
man observations by Falgarone et al. (2008) suggest a 
lower mass-to-flux ratio tor high-density regions, as we 
would expect. 
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3. RESULTS AND DISCUSSION 

3.1. Accretion History 

The accretion history for run E is shown in Figure [I] 
The first sink particle that forms accretes more than 
28 Mq during the simulation runtime of 0.683 Myr, even- 
tually becoming by far the most massive sink. Only after 
a delay of almost 20kyr after the formation of the first 
sink do secondary sinks begin to form. None of them 
gets more massive than 6M . In addition to the effect 
of radiation feedback, which raises the Jeans mass lo- 



cally and reduces fragmentation (see |Peters et ah] 2010b 
for a thorough discussion), fragmentation is further re- 
duced in this simulation by the magnetic field, just as 
is found in models of low-mass star formation (see Sec- 
tion 1). The same effects of reduced fragmentation occur 
even in our model of high-mass star formation with initial 
(M/$) » (M/$) cr . 

Figure [I] also shows the accretion rates and protostel- 
lar masses of the first sinks to form in run A, run B, 
run D, and run E for comparison. The accretion rate in 
the initial accretion phase is higher in run E than in the 
other simulations, but the long-term accretion behavior 
does not appear to be substantially different, apart from 
a stronger variation beyond 0.64 Myr. However, we note 
that in run B, fragmentation-induced starvation (Peters 
et al.||2010a|b ) terminates accretion onto a 23 M sink, 
whereas the 28 M sink in run E continues to accrete un- 
til the end of the simulation. The central star in run E 
can grow to larger masses because in the initial phase 
fragmentation is delayed by magnetic support. This al- 
lows the central object to maintain a high accretion rate 
for a longer time. Nevertheless, the accretion rate of 
the massive star in run E drops significantly when sec- 
ondary sink particles form, since they form in a dense 
ring around the central massive sink and to some extent 
starve it of material, even though they never cut off ac- 
cretion entirely. 

The stronger initial accretion phase in run E compared 
to run B yields the main contribution to the larger final 
mass of the massive sink. The magnetic field very effi- 
ciently redistributes angular momentum, resulting in an 
increased radial mass flux through the high-density equa- 
torial plane (see Section 3.2). Consequently, the initial 
accretion rate in the magnetized simulation (run E) even 
lies considerably above the non-magnetized single sink 
calculation (run A). 

The total accretion histories of all sink particles com- 
bined in each of the four simulations are contrasted in 
Figure [2] For the first 20kyr, the total accretion rate of 
run B, run D and run E is nearly identical. The accretion 
rate in these multiple sink simulations is generally higher 
than in the single sink run A since the large group of sinks 
accretes from a large volume, without needing to rely on 
outward angular momentum transport to deliver mate- 
rial to the direct feeding zone of the central sink. In the 
magnetized run E, no secondarysink particles form dur- 
ing the first 20kyr (see Figure IT]), so that the increased 
accretion rate in this phase is due to the additional an- 
gular momentum transport performed by the magnetic 
field. After the initial 20kyr, the total accretion rate in 
run E falls below run B and the control run D, deviating 
increasingly with time, but always staying above the ac- 
cretion rate in run A. The accretion histories of run B and 
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Fig. 1. — (a) Accretion history for run E. The first sink that 
forms eventually becomes the only massive one to form during 
the simulation runtime. Radiative heating and magnetic support 
reduce fragmentation and lead to the formation of only a small 
number of sinks in comparison to the equivalent gas dynamical 
model, (b) Accretion rates of the first sink particles in run A, 
run B, run D and run E. While the initial accretion phase for the 
first sink in run E leads to a higher accretion rate, the long-term 
evolution of the accretion rates looks very similar. In contrast to 
run B, however, in which accretion onto the first star stopped when 
it had reached 23 M after 0.70 Myr, accretion onto the 28 M 
star in run E still continues (compare (c)). The accretion rate 
drops significantly when secondary sink particles form in run E. (c) 
Protostellar masses of the first sink particles in run A, run B, run D 
and run E. Radiative heating (run B) and presence of magnetic 
fields (run E) increase the final masses of the massive stars. The 
largest part of the additional mass accretion in run E compared to 
run B is due to the stronger initial accretion phase. 
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Fig. 2. — Total accretion history of run A, run B, run D and 
run E. While the total accretion rates of run B and run D agree 
until ionization feedback stops massive star formation in run B, the 
total accretion rate in run E already starts to decline after 20kyr. 
The magnetic field in run E additionally supports the gas against 
collapse, reducing the total accretion rate. 

run D only start to separate at relatively late time when 
the ionizing radiation begins to terminate accretion onto 
the most massive s inks in run B (see the discussion in 
Peters et al.||2010b ), but the magnetic field additionally 
reduces the rate at which gas collapses in run E. Thus, 
the total accretion rate with magnetic fields is lower than 
wi thout magnetic fiel ds, in agreement with previous work 
by | Wang et al.| ( |2010| ) that used a stronger magnetic field. 

3.2. Magnetic Field Structure 

As the molecular cloud collapses, it forms a thin, ro- 
tationally flattened structure in the midplane of the ro- 
tating cloud. This disk-like structure is very similar to 
the one we find in the run without magnetic fields. As 
this structure does not have a Keplerian rotation curve, 
we avoid calling it a disk, though it behaves in many 
ways similar to one. This s tructure is es s entially what 
was called a pseudodisk by Galli fc Shul ( |1993a|b[ ). In 
the midplane, the dense, rotating gas drags the mag- 
netic field with it, winding the magnetic field lines into 
a tor oidal configuration that drives a magnetic tower 
flow dUchida fc Shiba^ 

sumoto fc Tomisaka 2004; Machida et al. 2004; Banerjee 
fc Pudritz|2006||Hennebelle fc Fromang|2008| ) However, 
the magnetic tower can only build up as long as the field 
lines are anchored in a coherently rotating flow. 

This velocity coherence is present initially but gets in- 
creasingly disrupted by fragmentation. Figure [3] shows 
slices of density, azimuthal velocity and vertical momen- 
tum through the midplane at three different times in 
run E. Initially, the whole region is rotating in the coun- 
terclockwise direction. As fragmentation proceeds, re- 
gions with clockwise (negative) azimuthal velocity ap- 
pear near the center of collapse, and these regions ex- 
pand radially with time. This means that the magnetic 
field lines that are anchored in the center cannot wind 
up anymore, and thus the vertical expansion of the mag- 
netic tower flow stalls. The plots of vertical momentum 
of the gas show that the magnetic outflow is thereafter 
mostly launched from a ring around the center whose ra- 
dius grows with time. The central regions, where the gas 



is rotating clockwise in some regions, eject little further 
material into the outflow. 

Figure [4] illustrates why gravitational fragmentation 
poses a substantial problem to maintaining the mag- 
netic tower flow. The figure shows density slices per- 
pendicular to the midplane through the magnetic tower 
together with the magnetic field. Shortly after gravita- 
tional collapse has formed a dense structure in the mid- 
plane and sink particles have formed, the field lines are 
nicely aligned and anchored in the rotating, dense gas, 
so that a magnetic tower starts to build up. However, as 
gravitational fragmentation proceeds, the circular mo- 
tion in the inner region gets increasingly disturbed by 
gravitational instability as well as angular momentum 
transport by the tower flow, resulting in a loss of veloc- 
ity coherence. Hence, the magnetic tower grows laterally 
into regions where the velocity field remains coherent as 
fast as it grows in height, so it is unable to produce a 
collimated outflow. This structure is a magnetic bubble 
rather than the outflow we are familiar with in low-mass 
star formation (I Cabrit & Andre 19 9l] Bachiller||1996 
|Reipurth fc BallypOOTj |Bally et al.pKM : 



Although the velocity field is not turbulent initially, 
a substantial amount of turbulence develops during the 
simulation runtime because of gravitational fragmenta- 
tion. Conversely, any initial turbulence wi ll die away in 
a sound-crossing time (M ac Low et al.|1998 ). We do sus- 
pect that initialization with highly supersonic turbulence 
might modify the large-scale dynamics and the shape of 
filaments in our flattened structure but not that it will 
qualitatively change our results. 

While it has been shown that magnetically driven out- 
flows from l ow-mass stars can still form in a turbulent en- 
vironment ( |Matsumoto fc Hanawa|[2010[ ), it is less clear 
whether highly collimated, magnetically- driven outflows 
will persist around massive stars since turbulent frag- 
mentation may even enhance the destructive effect of 
gravitational fragmentation discussed here on the veloc- 
ity structure of disks around high- mass stars. 

H II region formation further influences the magnetic 
field structure. Figure [5] shows slices of magnetic en- 
ergy density and the plasma beta (/3 P \ = Pth/Pmag with 
the thermal and magnetic pressure p t h and p m ag, respec- 
tively), along with magnetic field vectors. The left-hand 
plot demonstrates that the expansion of the H II region 
driven by ionization heating reduces the magnetic energy 
content inside the H II region of the most massive star 
in run E by more than four orders of magnitude. This is 
because the small amount of flux threading the initially 
ionized gas gets redistributed over the large volume that 
it expands into because of flux freezing. The alignment 
of the magnetic field, which remains visible outside the 
H II region, also gets destroyed. The right-hand plot 
shows the same region after the former H II region has 
completely recombined. The magnetic field energy has 
only increased slightly, since the gas has not been com- 
pressed much, and the magnetic field orientation is still 
disarranged. 

3.3. Gravitational and Magnetic Torques 

To better understand the reasons for the breakdown of 
coherent rotation and the onset of counter-rotation in the 
flattened accretion flow, which ultimately disrupts the 
magnetic outflow, we study the different torques acting 
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Fig. 3. — Slices of gas density p, (left) and azimuthal velocity (middle) in the midplane and (signed) vertical momentum density pv z 
(right) slightly above the midplane for three different snapshots of run E. The mass of the most massive star in the cluster is given in the 
images. Fragmentation sets in after the formation of the first filaments and sink particles (black dots). The resulting torques can reverse 
the direction of the azimuthal velocity component. The arrows indicate the direction of velocity. Although the cloud is initially rotating 
counterclockwise (white), the central region soon begins to rotate clockwise (red), and this region of reversed azimuthal velocity gradually 
expands. This has consequences for the driving of the magnetic outflow, visible in the slices of vertical momentum density (right). The 
outflow is primarily launched from its outer boundary, where coherent rotation is maintained, while momentum transport from the inner 
region is relatively weak. The online material contains an animated version of this figure. 
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Fig. 4. — Magnetic bubble and field structure. The images show density slices parallel to the rotation axis and perpendicular to the 
midplane through the magnetic bubble. The arrows indicate the direction of the magnetic field lines (top) and velocity (bottom), and the 
black dots show the positions of sink particles. In the initial phase (left), the field lines wind up to build up a magnetic tower flow. As 
the fragmentation occurs, the velocity coherence in the central regions disappears, so that magnetic field lines anchored there stop winding 
up (right). Hence, the magnetic bubble grows primarily laterally but does not form a more collimated outflow. Note the different spatial 
scales in the images. 
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Fig. 5. — Influence of ionizing radiation on magnetic field structure. The images show slices of magnetic energy density e m ag (top) and 
the plasma beta f3 p \ = Pth/:Pmag (bottom) along with magnetic field vectors and sink particles. The ionizing radiation emitted by the most 
massive star of the cluster creates an H II region, which dramatically decreases the magnetic energy content in this region and destroys 
the coherence of the magnetic field structure (left). After complete recombination of this H II region, the magnetic energy density remains 
reduced and the magnetic field structure is not reconstructed (right). The black dots represent sink particles. The online material contains 
an animated version of this figure. 



on the gas in the central regio ns. The torque r exerted on 
the gas in the midplane (e.g., Ba nerjee fc Pudritz| |2006 ) 
is given by 



tional torques 



-"J 



' mag 



(3) 



J r x [(V x B) x B] dV 



(4) 



with the spherical radius r, the radial velocity v r , the 
cylindrical radius R, the azimuthal velocity and the 
differential solid angle d£l. Since the purely hydrody- 
namical simulation (run B) and the MHD calculation 
(run E) show fragmentation in the midplane, the nega- 
tive torques required to reverse the direction of rotation 
should be present in both cases. Figure [6] demonstrates 
that his is indeed the case. The braking and reversal of 
the coherent velocity field is thus not a new phenomenon, 
but there are some notable differences in the fragmenta- 
tion behavior. 

While the only torques acting in the hydrodynamical 
case are gravitational, the magnetic field can exert addi- 



on the gas in the MHD case. Under the assumption of 
cylindrical symmetry, the torque exerted by the magnetic 
field is 

1 



' mag 



47T 



BrRBp dtt 



(5) 



with the radial magnetic field component B r and the 
toroidal component B^. The difference 

(6) 



' grav 



' mag 



can be attributed to gravitational torques. In principle, 
both r mag and r grav can become negative and force the 
velocity field to reverse. 

Since fragmentation breaks cylindrical symmetry, it is 
questionable whether Equation ([5| can still be applied 
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Fig. 6. — (Left) The cumulative torque r is plotted for the hydrodynamical simulation (run B) and for the MHD simulation (run E). 
Negative torques, indicating braking and reversal of the azimuthal velocity field, appear in both cases. (Right) Cumulative magnetic 
torques calculated in two different ways, the z-component of Equation Q and via Equation |5j. Deviations become visible only for radii 
> 10 17 cm, when the averaging volumes become disparate. The snapshot is taken from the end of run E, where the rotationally flattened 
flow is strongly fragmented and asymmetric. 



to calculate the magnetic torque. To test this, we have 
compared r mag computed via Equation ([^ with r mag ^, 
the ^-component of r mag from Equation Ml). In the lat- 
ter case, we have evaluated r mag ^ on a Homogeneous 
grid with a 98 AU cell size, corresponding to the high- 
est refinement level in the adaptive mesh, of dimensions 
0.49 x 0.49 x 0.06 pc 3 . A typical snapshot is shown in Fig- 
ure [6] The agreement is surprisingly good. Although the 
data is not strictly cylindrically symmetric, the curves 
agree very well after the angular averaging process, at 
least for radii smaller than ~ 10 1T cm, where the regions 
over which the data is averaged in the two calculations 
differ little. Hence, Equation ([5| is still applicable to 
study the radial dependence of r mag . 

The shape of the filaments is substantially different 
in run B and run E. While the filaments are randomly 
oriented in the purely hydrodynamical case (run B), the 
magnetic field organizes the filaments into ring-like struc- 
tures in run E, as can be nicely seen in the middle panels 
of Figure [3] This raises the question whether the neg- 
ative magnetic torques dominate over the gravitational 
torques in the MHD simulation. This is generally not 
the case, as is shown in Figure [7| Although there are re- 
gions where magnetic and gravitational torques are close 
to equipartition, the gravitational torques are generally 
larger in magnitude than the magnetic ones. Hence, it is 
commonly the gravitational torques that are responsible 
for braking and velocity reversals, in agreement with the 
observation that we find these effects also without any 
magnetic fields in run B. Magnetic braking appears to be 
globally less relevant compared to gravitational torques. 
However, for the dynamics in the inner region and the ac- 
cretion onto individual sink particles from a smaller-scale 
accretion disk, local magnetic braking can nevertheless 
be important. 

Figure [7] also shows that the negative total torques can 
be associated with circularly arranged filaments in the 
midplane (marked as circles). Since the filaments often 
represent the boundary between clockwise and counter- 



clockwise rotation (compare Figure [3k, this is exactly 
where negative torques are expected. When the filaments 
are circularly arranged, the angular averaging can yield 
pronounced negative peaks in the radial torque distribu- 
tion. Of course, this special condition is not the only one 
that can lead to negative torques, but it regularly occurs. 

The velocity vectors plotted in Figure [7] indicate that 
the negative torques cause a reversal of the radial ve- 
locity component at the filaments. This reversal of the 
radial velocity component is not simply the effect of an 
expanding motion near the center, but a result of the 
gravitational torques exerted by the non-uniform gravi- 
tational instability in the disk. 

3.4. Gravitational Instability and Fragmentation 

The fragmentation of the rotationally flattened accre- 
tion flow leads to the formation of relatively long-lived, 
circularly arranged filaments in run E. This is in contrast 
to the situation in run B, where the filaments are gen- 
erally less elongated and their shape and location can- 
not directly be associated with the rotating flow (see 
Figure [8]). To understand the apparent stability of this 
structure, we investigate the relative importance of two 
contribut ing m echanisms: support by magnetic pressure 



(Section 3.4.1) a nd sup port by strong shear across the 



filaments (Section 3.4.2 ). 



3.4.1. Support by Magnetic Pressure 

Though the magnetic field appears to be not the main 
driver of negative torques, it is clearly dynamically rel- 
evant for the gas dynamics within the accretion flow. 
This is demonstrated by examining the plasma beta 
Pp\ = Pth/.Pmag in the midplane, as shown in Figure JSJ 
Outside of the H II region around the massive star, trie 
magnetic pressure p mag exceeds the thermal pressure p t h 
generally, except in the dense filaments. In the central 
region, p mag is more than two orders of magnitude larger 
than p t h, but near the filaments p mag and p t h reach equal- 
ity. 
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Fig. 7. — Torques and density structure. The figure shows cumulative total, magnetic and gravitational torques as function of radius 
and the gas density p in the midplane for two representative snapshots. The magnetic r m ag and gravitational torques r gr av can reach 
equipartition locally, but generally the gravitational torque is considerably larger in magnitude. Regions of negative total torques (marked 
as lines) can be associated with circularly arranged filaments (marked as circles) at the same radii. The black dots represent sink particles. 
The arrows indicate the direction of the velocity field. 



The plasma beta f3 p \ also quantifies the importance 
of magnetic support against gravitational collapse in a 
linear perturbation analysis. The Toomre Q-parameter 
is defined as 

« = Js < 7 > 

with epicyclic frequency ft, sound speed c s , surface den- 
sity E, and Newton's gravitational constant G. Linear 
gravi t ational instability sets in when Q < 1 (| Toomre 



1964 Goldreich & Lynden-Bell 1965). When mag netic 



fields are taken into account ( |Kim & Ustriker|2001| ), one 
can define the magnetic Toomre parameter 



3.4.2. Support by Shearing Motion 

Since the Toomre analysis as presented above only ap- 
plies to linear perturbations to Keplerian rotation, it en- 
tirely neglects the strong shear along the filaments (com- 
pare Figure |3| . Because the gas on one side of a filament 
is moving clockwise while the gas on the other side is 
moving counter-clockwise, there is a large velocity gra- 
dient across the filament. This velocity gradient appears 
able to stabilize the filament. To see how big this effect 
is, we define the shear timescale 



1 



Qm = 



K (c 



a/ 2 



Vxvl 



(10) 



ttEG 



(8) 



with the vorticity V x v and compare this with the free- 
fall time scale 



with the Alfven velocity va- Since /3 p i 
of Q and Qm is 



- s2 




c sl v A-> ^ ne ratio 



(9) 



iff 



3tt 
Z2Gp' 



(11) 



Because (3 P \ is larger than unity in the filaments, Qm and 
Q only differ by a factor of order unity there, so that 
magnetic pressure plays no crucial role in stabilizing the 
filaments. 



The ratio t^/t^ is plotted in Figure [8] The filaments 
lie in the marginally stable regime, so that this simple 
timescale comparison is not totally conclusive. There is, 
however, a clear difference between the organization of 
the filaments in run B and run E visible in Figure [8] 
In run E, the filaments get stretched by the shear and 
organize circularly, while the filaments in run B are much 
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Fig. 8. — (Left) Slice of gas density in the hydrodynamical simulation run B. The filaments are much more roundish and are not stretched 
between rotating and counter-rotating gas. The black dots represent sink particles. The contours show the density structure. (Middle) The 
plasma beta f3 p \ = i?th/Pmag in the midplane of run E. Outside the H II region around the 26 Mq star, the magnetic pressure is generally 
larger than the thermal pressure and hence is dynamically relevant. However, p t h is slightly larger than p m &g in the filaments, and thus 
magnetic pressure plays no dominant role in stabilizing them. (Right) The ratio £ff/£ s h f° r the same region from run E. The filaments, 
marked with the density contours, are on the boundary between the stable and the unstable regime. 



more roundish in shape. In run E, the filaments separate 
rotating and counter-rotating gas, which is not generally 
the case in run B. Since it is the magnetic field that 
organizes the filaments in this way, the magnetic field 
indirectly stabilizes the filaments, though not directly 
through magnetic pressure. 

3.5. Relevance for Protostellar Jets 

These findings are immediately applicable to the large- 
scale rotating, infalling flow, far away from the ac- 
tual radius where high- velocity protostellar jets must be 
launched. Since the rotation of the launching region de- 
termines the outflow velocity, our larger-scale magnetic 
tower flow, with peak velocities around 5 km s _1 , (com- 
pare Figure [4| is slower by a factor of more than 20 
than optical jets w ith velocities exceeding 100 km s _1 (eg. 
|Bally et al.|[2QQ7| ). 

However, the behavior seen here may have important 
implications for understanding the role of jets in high- 
mass star formation. Such jets are driven from radii of 
less than an astronomical unit in the accretion disks of 
low mass stars. We do not resolve the accretion flow to 
these small radii in our models. However, the huge accre- 
tion rates required for high mass star formation suggests 
that those inner regions must also be highly perturbed 
and may likely be gravitationally unstable, and thus sus- 
ceptible to the same fragmentation and destruction of 
velocity coherence that we have seen in our mo dels. Ob- 
servations of a counter-rotating accretion disk (Remijan 



such small scales if the magnetic field outside this radius 
is totally disarranged. Models of dynamo generation of 
fields in such disks including non-ideal MHD effects are 
required to understand whether coherent fields capable 
of driving jets can form on such small scales. Our find- 
ings thus raise the question whether highly collimated, 
magnetically-driven jets from massive protostars can sur- 
vive. A negative answer could also explain why out flows 
from massive stars appear to be less collimated ( |Arce| 



2005): they are driven 



Hollis 2006) demonstrate that this pos sibility ex ists. 
Simplified one-dimensional calculations (Vaidya et al. 



2009) indicate that the temperature may grow strongly 



enough below a radius of 1 AU to stabilize the disk, but 
it may be difficult for the magnetic tower to build up on 



et al ||2007| |Beuther fc Shepherd! 

by the ionizing radiatio n instead ol cohere ntly rotating 
magnetic fields. Indeed, Peters et al. (2010a) presented 
a detailed comparison between the ionization-driven out- 
flows in the purely hydrodynamical simulation with ob- 
servations of W51e2 and found that the molecular and 
ionized parts of the outflow have a very similar velocity 
structure, both qualitatively in the characteristic features 
as well as quantitatively in the magnitude of the velocity 
gradient. 

3.6. Morphology of H II Regions 

To study the morphological types of the H II regions in 
run E in more detail, we have generated synthetic VLA 
observ ations at 2 cm, follow ing the algorithm described 
earlier ( Peters et al.|2010a|c| . Figure [9jshows H II regions 
of all morphological types round in run E. There seems 
to be no systematic difference between the morphologies 
in the magnetic and purely hydrodynamic simulations 
( |Peters et ; aLl2010a|c| >. 

Although trie magnetic field does not appear to have 
any influence on the H II region morphologies themselves, 
it does play a significant role in determining their size. 
Figure [10] shows a comparison of thermal and magnetic 
pressurefor the H II region from Figure [5] The thermal 
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Fig. 9. — H il region morpholgies in run E. The figure shows synthetic maps of free- free emission from ultracompact H II regions around 
the massive protostar at different time steps and from different viewing angles. The cluster is assumed to be 2.65 kpc away, the full width 
at half maximum of the beam is / /14 and the noise level is 10 — 3 Jy. This corresponds to typical VLA parameters at a wavelength of 2 cm. 
The protostellar mass of the central star which powers the H II region is given in the images. Black dots represent sink particles. 



pressure inside the H II region is of the same order as 
the magnetic pressure immediately outside the H II re- 
gion. Thus, the magnetic pressure yields an important 
contribution to the total pressure in the environment of 
the H II region, constraining its expansion. We therefore 
expect H II regions in the presence of strong magnetic 
fields to be generally smaller than without the magnetic 
field. 

4. SUMMARY AND CONCLUSION 

We have presented the first three-dimensional, 
radiation-magnetohydrodynamical collapse simulation of 
massive star formation including heating by both ioniz- 
ing and non-ionizing radiation. We used sink particles 
to represent accreting protostars. We compared the ac- 
tion of the magnetic field in this high-mass star forma- 
tion simulation to the major ways that it acts in low- 
mass star formation. As we summarized in the introduc- 
tion, these include giving support against gravitational 
collapse, transferring angular momentum outward, and 
driving outflows. 

We find that, although we started with a mass-to-flux 
ratio that was more than an order of magnitude super- 
critical, (M/$) = 14(M/$) cr , the magnetic field still can 
reduce fragmentation (Section 3.1 ), just as it does in low- 
mass star formation. The additional magnetic support 
prevents the gas from collapsing into as many secondary 
fragments compared to the non-magnetic case, leaving 
the most massive protostar with a larger gas reservoir. 
For the same reason, the total accretion rate is reduced. 
However, it is important to note that the magnetic field 



does not suppress fragmentation completely. The num- 
ber of fragments is reduced at most by a factor of two. 
Similar findings are reported by Hennebelle et al. (2010, 
A&A, submitted). 

Additionally, magnetic braking near the center reduces 
the angular momentum of the collapsing gas even further. 
Thus, more material is transported radially inwards onto 
the massive central protostar. This is important since the 
distribution of mass accretion among the protostars in 
the fragment at ion- induced starvation scenario depends 
sensitively on the (in) efficiency of gas transport tow ards 
the center by internal torques ( |Peters et al.||2010a|b[ ). 

The combination of these two effects causes the for- 
mation of a 50% more massive central star compared to 
our equivalent simulation without magnetic fields (Fig- 
ure [Tp) . If th e magnetic field strengths inferred by Fal- 
garone et al. (|2QQ8[) are typical throughout high-mass 



clumps, and not just in their high-density cores, mag- 
netic braking and reduction of fragmentation may be 
even stronger than reported here. 

The winding up of the magnetic field into a toroidal 
configuration leads to the formation of a large-scale tower 
flow, surrounding all the protostars. However, two effects 
tend to weaken and broaden the outflow. The first one 
is the gravitational fragmentation of the accretion flow. 
This disrupts the circular motion necessary to drive the 
tower flow in the central region. Second, the thermal 
pressure of ionized gas is by far larger than the magnetic 
pressure and hence dynamically dominant within the H II 
region. The resulting expansion of the H II region dra- 
matically reduces the magnetic energy content within it, 
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Fig. 10. — Comparison of thermal and magnetic pressure for the data from the left-hand panels in Figure [5] The thermal pressure p t h 
inside the H II region (left) is of comparable magnitude to the magnetic pressure p ma g outside the H II region (right). Thus, magnetic 
pressure plays a significant role in constraining the size of expanding H II regions. The black dots represent sink particles. 



and tangles the field lines, further weakening the flow. 

The torques required to brake the large-scale rotat- 
ing flow and reverse the azimuthal velocity are pri- 
marily gravitational, though magnetic and gravitational 
torques can reach equipartition locally. Magnetic brak- 
ing appears to be of minor importance for the global 
gas dynamics, but not necessarily for the local accretion 
onto protostars. Again, if stronger magnetic fields are 
imposed, magnetic braking might become qualitatively 
more important. The shape and dynamics of the fila- 
ments is vastly different in the magnetohydrodynamical 
(run E) and purely gas dynamical (run B) case: while 
the filaments in run B are disorganized, roundish and 
quickly collapse to form stars, the filaments in run E are 
stretched between clockwise and counter-clockwise rotat- 
ing gas and can be maintained over a large fraction of the 
simulation runtime. 

If the fragmentation process observed in our simula- 
tions holds down to the smallest scales of the massive ac- 
cretion disk, as suggested by the very high accretion rates 
of the order lO -3 M0yr -1 inevitably needed to form a 
massive star, the magnetic field can only wind up at radii 
less than 1 AU, where the temperature is high enough 
to suppres s gravitationl instability of the disk ( |Vaidya| 
et al. 2009). In this case, magnetically-driven, steady jets 
around massive protostars can only survive until gravita- 
tional fragmentation disrupts uniform rotation, and ion- 
izing radiation becomes dynamically relevant. We call 
this process fragmentation-induced outflow disruption. 
The fast jets launched from the inner disk region should 
then be highly episodic like the accretion rates. The 



uncollimated outflows from massive stars may be better 
explained as driven by the ionization feedback. 
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